Comment on "Regularizing capacity of metabolic networks" 
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In a recent paper, Marr, Muller-Linow and Hiitt [Phys. Rev. E 75, 041917 (2007)] investigate an artificial dynamic 
system on metabolic networks. They find a less complex time evolution of this dynamic system in real networks, 
compared to networks of null models. The authors argue that this suggests that metabolic network structure is a 
major factor behind the stability of biochemical steady states. We reanalyze the same kind of data using a dynamic 
system modeling actual reaction kinetics. The conclusions about stability, from our analysis, are inconsistent with 
those of Marr et ol. We argue that this issue calls for a more detailed type of modeling. 
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[. INTRODUCTION 

Within living organisms, matter is constantly converted be- 
tween different molecular species. It is often assumed that 
concentrations of metabolites tend to settle into steady states 
(rather than showing periodic or chaotic behavior) (Q]). Fur- 
thermore, experimental studies of metabolic pathways are typ- 
ically performed under steady-state conditions O). Although 
the steady state is, strictly speaking, a mathematical abstrac- 
tion, it is nevertheless a useful reference state (3). The steady- 
state assumption is also fundamental to traditional metabolic 
control analysis (2). In a recent article (4) Marr, Muller-Linow 
and Hiitt hypothesize that the network structure of metabolism 
is an important factor for promoting a steady state, rather 
than a complex, dynamics. The authors run an artificial dy- 
namic system on the networks. In short, a vertex i can have 
two states, or 1. If the sum of the state variables in the 
neighborhood of i exceeds a fixed threshold, then i changes to 
the other state. This scheme is very different from real bio- 
chemical dynamics and traditional models of reaction kinet- 
ics. First, the state variables in metabolic models are usually 
continuous (concentrations (0)) or sometimes discrete vari- 
ables (molecule counts (01)), but never binary. Second, the 
sum of these variables is conserved (if in- and outflow is ne- 
glected); whereas, in the dynamics of Marr et al., this is not 
the case. This dynamics does in fact not reach a steady state; 
instead, the authors analyze the complexity of the output time 
series with entropy-like measures. They conclude that real 
metabolic networks give a less complex output than different 
ensembles of model networks, and argue that this implies that 
the structure of real metabolic networks may promote steady 
state dynamics. However, no theory for how the values of the 
"entropies" of the binary dynamics relate to the stability of 
metabolic flux is given. Many authors have used binary dy- 
namics to model other processes of cellular biology, such as 
signal transduction (e.g. Ref. (6)) or genetic regulation (e.g. 
Ref. (|7|)). However, these models explicitly try to describe 
systems in which the components, at least under some circum- 
stances, show switch-like, binary behavior (6). Moreover, the 
lack of detailed understanding of e.g. gene expression makes 
it more natural to use simplified dynamics to uncover under- 
lying principles. Metabolic reactions follow known principles 
of chemistry and are described by well-established differential 




FIG. 1 (Color online.) Illustration of the construction of substrate 
graphs from chemical reactions (a), and our method for recreating a 
reaction system from a substrate graph (b). 



equation models. Therefore, in contrast to the mentioned in- 
formation processes, there is no need for, or natural interpre- 
tation of, a binary description of the metabolism. Therefore 
the conclusions of Marr et al. can be validated by comparing 
their results to results obtained using such standard differen- 
tial equation modeling. We follow this approach to find no 
significant difference between real and null-model network. 
Our simulations in this Comment are rather limited, but suf- 
ficient to make the point that Marr et al. 's study cannot rule 
out the possibility that stable steady-states can be explained 
more simply, as direct consequences of the reaction kinetics, 
rather than as a result from the interaction between the net- 
work topology and the dynamic system. 



II. CONSTRUCTING A REACTION SYSTEM 

Marr et al. 's examples of real metabolic networks are taken 
from a work of Ma and Zeng (MZ) (|8|). They are substrate 
graphs where a substrate is linked to the products of a re- 
action (Fig[T}a)). MZ also preprocessed the data by omitting 
ubiquitous "currency metabolites" (10). We also use the MZ 
networks as the starting point for our simulations. In con- 
structing the networks, information is lost; so if one wants to 
simulate the reaction system with a substrate graph as start- 
ing point, one needs (explicitly or implicitly) to recreate the 
reactions via a model. From two assumptions about reaction 
systems, we propose a simple scheme to create a plausible set 
of reactions that can be reduced to a given substrate graph. 
We first assume that all reactions are of a 2-2-form: A + B <-> 
C + D, or 2-1 -form: A + B «-» C. These are the most com- 
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mon forms of biochemical reactions. We do not include more 
complex reactions, mostly because it would significantly in- 
crease the computational complexity. Our second assumption 
is that the number of reactions creating an edge in the sub- 
strate graph is rather small (we confirm a posteriori that the 
average number of reactions per vertex is similar to that of the 
real substrate graph). The following algorithm is a simple way 
of fulfilling these assumptions: 

1 . Start with all edges unmarked. 

2. Pick one unmarked edge (A, C). 

3. Find the maximal, full bipartite subgraph K(A,c) (a sub- 
graph consisting of two vertex sets, and edges between 
every pair of vertices in the different sets, but no edge 
between vertices in the same set) that contains (A, C). 
(See Fig. Eb).) 

4. Pick one four-cycle of ^T(a,c) including (A, C) (say 
(A, C, B, D, A)). (See Fig. [lib).) 

5. Add A + B<-»C + Dto the set of reactions and mark 
the edges (A, C), (C, B), (B, D), (D, A). 

6. If there are unmarked edges go to step [2] 

K(a,C) gives all reactions of the 2-2-form that induce the edge 
(A, C) in a substrate graph. If ^T(a,c) 1S empty at step [3] then 
(A, C) must have been generated by a reaction of a 2-1 -form. 
Therefore, instead of #(a,c)> consider the set L(a,c)> of (A, C) 
and its adjacent edges, and deduce reactions of the 2-1 -form 
at step [5] 



III. REACTION KINETICS 

Once we have a set of reactions, derived from a sub- 
strate graph, we simulate the biochemical dynamics by stan- 
dard Michaelis-Menten (MM) kinetics (fill) supplemented by 
noise in the enzyme and substance concentrations. MM ki- 
netics are used to model enzyme-catalyzed reactions, which 
are common in metabolism. The MM description builds on 
the principle of mass-action — the forward rate of a reac- 
tion is proportional to the product of the concentrations of the 
substrates — and uses additional assumptions to simplify the 
resulting rate laws. For the reaction i, A + B — » C + D (which 
is assumed to be enzyme-catalyzed), we use a version of the 
two-substrate MM rate law to calculate the flux 
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where t/e is a random variable modeling fluctuations in the en- 
zyme concentration. The time evolution of the concentration 
of a substance A is then determined by 
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FIG. 2 (Color online.) Time evolution of the average standard de- 
viation £ of the flux. Standard errors are smaller than the symbol 
size. 



77 conc . is a random noise term modeling fluctuations due to in- 
and outflow of A. We seek to use the same information as 
in Ref. (4), therefore we do not use empirical parameter val- 
ues (which are, anyway, hard to obtain). Instead, we assign 
parameters in arbitrary units, but we choose them to be rea- 
sonable relative to one another. t]e and ?7 conc . are normally 
distributed A^(0, 0.002) (the first argument is the mean; the 
second is the standard deviation). New values for t]e and 
i]conc. are drawn every time step of the integration. The initial 
values of substance concentrations, enzyme concentrations 
and reaction coefficients (ky, k\, k 2 and k\ 2 ) are drawn from 
N(l, 1), iV(0.2,0.2) and N(0.1,0.1) respectively. We choose 
fi k = 1 x 10~ 3 , (r k = 5 X 10~ 4 , fi = 1 and cr = 0.5. The 
equations are integrated with a second order Runge-Kutta- 
Helfand-Greenside scheme with time step 0.1, total running 
time 2500 time units, and 20 averages over different sets of 
initial configurations. These runs are, for each network, aver- 
aged over 20 realizations of the reaction-system construction. 
For comparison, we also run the dynamics on 100 samples of 
one Marr et al. 's null-model networks — random graphs with 
the same degree sequence as the original network. These are 
obtained by random rewiring — we go through all edges and 
for each edge (i, j) randomly pick an edge (/', /) and replace 
these edges with (i, /) and (i', j) (unless this would introduce 
a multiple edge or a self-edge, in which case a new random 
edge (f, f) is selected). 



IV. SIMULATION RESULTS 

In the simulations, a vast majority of the substance con- 
centrations converge to steady states. A few subnetworks of 
oscillating or chaotic concentrations may exist but in this 
Comment we focus on bulk properties. To study the approach 
to equilibrium, and fluctuations, we measure the average stan- 
dard deviation of the flux through the substances, 



(2) 



where 



(3a) 
(3b) 



where the sum is over all A's reactions, the sign in front of 
r,- is positive (negative) if A is a product (substrate) of i, and 



The sums in Eq.[3a]are over all substances; the sum in Eq.[3bl 
is over fs reactions; the factor 1/2 comes from the double 
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count of mass-flow in and out of substance i, and n is the num- 
ber of substances. In a plot of £(f), the stability of the steady 
state can be monitored by two quantities. First, when the sys- 
tem approaches an equilibrium, "Lit) will decrease — a faster 
decrease implies a more stable system. Second, a lower equi- 
librium level means that the system has less cyclic or chaotic 
components, and responds faster to perturbations from the 
noise, and therefore has more stable steady states. In Fig. [2] 
we plot S(f) for MZ's human metabolic network and the null- 
model networks. The two curves almost overlap. For different 
null-model realizations the curves may deviate slightly from 
those of the real network, but there the null model cannot be 
rejected with any high level of significance. The equilibrium 
level is of the order of the noise, which means that periodic 
and chaotic behavior is almost fully suppressed. To be more 
systematic, we note that 2(f) display two aspects of stability — 
how fast equilibrium is reached and the height of the equilib- 
rium level. Since the curves do not fit any simple functional 
form we measure the half-time (the time to reach midway be- 
tween X(0) and X(oo)) using spline interpolation. The p-value 
(fraction of null-model observations lower than the real value) 
of the half-time is 67%; the corresponding value of the equi- 
librium level is 65%. The effect of the difference in network 
structure between the real and null-model networks is thus, 
in this case, negligible. We test a few other sets of parame- 
ter values, organisms (E. coli and M. musculus), and a more 
straightforward mass action kinetics; and, in all cases, arrive 
at the above conclusion. A full scan of the parameter space 
would be interesting, however, in this Comment we just make 
the point that the conclusion of Marr et al. can be inconsistent 
with more realistic simulations, and leave the insignificance 
of network structure to the stability of metabolic steady states 
as a conjecture. 



V. CONCLUSIONS 

To conclude; simple, stylized dynamic systems — "dynamic 
probes" — are valuable tools for studying complex biological 
systems. We believe that these should be designed to model 
the real dynamics as closely as possible. Marr et al. (3) pro- 
pose a dynamic probe to study metabolic steady states that 
violates many of the known features of reaction kinetics; it 
does not even reach a steady state — the objective of the study 
as given in their Abstract. By the collective effort of re- 
searchers, our understanding of the metabolism continuously 
advances. This does, however, not improve the approach of 
Marr et al. as their dynamics does not make use of biochemi- 
cal information. Relative to biological information processes, 
metabolism is a rather simple system, and is believed to be 
well-described by simple differential equation models. For 
the simulations we carry out, the null-model networks are not 
less efficient than the real networks in suppressing complex 
dynamics. If this holds in general, then steady-state stability 
is a fundamental property of chemical reaction systems. The 
study of Marr et al. does not rule out such a simple explanation 
of steady-state behavior — an explanation that is a common 
opinion in biochemical literature (cf. Hofmeyr and Cornish- 



Bowden's dictum "mass-action is the intrinsic driving force 
for self-organization of reaction networks" yl)). If the reac- 
tion kinetics is the sole cause of metabolic steady-states, the 
steady-state dynamics is a constraint to, rather than an out- 
come of, natural selection. This situation is reminiscent of the 
power-law degree distribution of metabolic networks. Such 
distributions can also be seen in astrochemical networks that 
are not subject to natural selection ( 13). We believe the ques- 
tion of dynamic stability in metabolism should be studied at a 
more detailed level than networks, avoiding the reduction of 
database information to substrate graphs. This is not to say 
that graph theory is useless in the study of metabolism. On 
a large scale, metabolic networks are different from the null- 
model networks. This is a valid conclusion from Ref. (0) — the 
authors manage to separate the real metabolic networks from 
both the null-model network we use, and various other types 
of model networks derived from the original graph. We be- 
lieve much information about the organization of metabolism 
lies in the answer to how this separation occurs. On smaller 
scales, network theory can be used to find e.g. functional mod- 
ules ( 10) and functions of individual metabolites dl4l) . 
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